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1.0 INTRODUCTION 


NASA’s advanced engine programs are aimed at progressively higher efficiencies, greater 
reliability, and longer life. Recent studies have indicated that significant engine performance 
advantages can be achieved by employing advanced seals [1]*, and dramatic life extensions 
can also be achieved. Advanced seals are not only required to control leakage, but are 
necessary to control lubricant and coolant flow, prevent entrance of contamination, inhibit the 
mixture of incompatible fluids, and assist in the control of rotor response. 

Recognizing the importance and need of advanced seals, NASA, in 1990, embarked on a five- 
year program (Contract NAS3-25644) to provide the U.S. aerospace industry with computer codes 
that would facilitate configuration selection and the design and application of advanced seals. 

The program included four principal activities: 

1. Development of a scientific code called SCISEAL, which is a Computational Fluid 
Dynamics (CFD) code capable of producing full three-dimensional flow field 
information for a variety of cylindrical configurations. The code is used to enhance 
understanding of flow phenomena and mechanisms, to predict performance of 
complex situations, and to furnish accuracy standards for the industrial codes. The 
SCISEAL code also has the unique capability to produce stiffness and damping 
coefficients that are necessary for rotordynamic computations. 

2. Generation of industrial codes for expeditious analysis, design, and optimization of 
turbomachinery seals. The industrial codes consist of a series of separate stand-alone 
codes that were integrated by a Knowledge-Based System (KBS). 

3. Production of a KBS that couples the industrial codes with a user friendly Graphical 
User Interface (GUI) that can in the future be integrated with an expert system to 
assist in seal selection and data interpretation and provide design guidance. 

4. Technology transfer via four multiday workshops at NASA facilities where the results 
of the program were presented and information exchanged among suppliers and users 
of advanced seals. A Peer Panel also met at the workshops to provide guidance and 
suggestions to the program. 

This final report has been divided into separate volumes, as follows: 

Volume 1: Executive Summary and Description of Knowledge-Based System 

Volume 2: Description of Gas Seal Codes GCYLT and GFACE 

Volume 3: Description of Spiral-Groove Codes SPIRALG and SPIRALI 

Volume 4: Description of Incompressible Seal Codes ICYL and IFACE 

Volume 5: Description of Seal Dynamics Code DYSEAL and Labyrinth Seal Code KTK 

Volume 6: Description of Scientific CFD Code SCISEAL. 

This volume describes the codes ICYL (incompressible, cylindrical) and IFACE 
(incompressible, face) that determine performance of a variety of incompressible fluid-film 


Numbers in brackets designate references presented in Section 6.0. 
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seal configurations. Both were written in FORTRAN for a PC environment using OS/2 as an 
operating system. References 2 and 3 provide the details of code implementation. 

The computer code ICYL was developed [1] to evaluate the performance of cylindrical seals 
operating with incompressible fluids. The computer code IFACE was subsequently developed 
to also handle face seals. The pressures generated in plain cylindrical seals with 
incompressible fluids typically result in forces which are normal to the displacement and 
therefore tend to destabilize the rotating shaft. Surface roughness, geometry alterations, and 
external pressurization are ways in which the direct stiffness and damping coefficients can be 
improved and the cross-coupled stiffness decreased in order to improve stability. 

The pressure and velocity distributions within the seal clearance are first evaluated from the 
governing equations. From these, design quantities such as seal leakage flows, power loss 
and resulting forces and moments are calculated. Minimum film thicknesses and maximum 
pressures as well as critical rotor-dynamics coefficients such as stiffness, damping and critical 
mass are evaluated. 

Program capabilities: 

1. 2-D incompressible isoviscous flow in cylindrical and face geometry. 

2. Rotation of both rotor and housing. 

3. Roughness of both rotor and housing. 

4. Arbitrary film thickness distribution, including features such as steps, pockets, tapers 
and preloaded arcs 

5. Rotor position described by four degrees of freedom (translational and rotational) 

6. Up to 32 dynamic coefficients as well as the critical mass may be calculated for use 
in rotor-dynamic design, including system response and stability calculations. 

7. Steady external forces and/or moments may be prescribed to find the position of the 
rotor relative to the housing. 

8. Pocket pressures or orifice size are prescribed. 

9. Laminar or turbulent flow. 

10. Cavitation. 

11. Inertia pressure drop at inlets to fluid film (from ends of seal and from pressurized 
pockets). 

Assumptions: 

1. The film thickness is assumed to be small compared with seal lengths and diameters 
but large compared with surface roughness. 

2. Pockets supplied from an external pressure source through an orifice restriction are 
assumed to be sufficiently deep that the pressure is constant within them. 

3. Wall roughness is assumed to be isotropic and represented by an "equivalent sand 
roughness" height. 

4. Fluid inertia effects in the film are negligible. 

This report describes the fundamental theory, sample problems and code validation. 
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2.0 THEORETICAL DESCRIPTION AND NUMERICAL 
METHODS 


Figures 1 and 2 illustrate the geometry of a cylindrical seals as well as the coordinate system 
used to describe it. Figure 1 shows the seal housing of length L separated from the rotor by 
the Film thickness C. The coordinate system is placed at the mid-length of the seal with the 
circumferential coordinate 0 measured from the x-axis. Figure 2 shows an axial cross-section 
of the film thickness with an eccentric rotor. Figure 3 similarly illustrates the geometry and 
coordinate system used to describe a face seal, from side and end views. The top part 
illustrates the seal ring rotating with a stationary mating ring, while the bottom illustrates the 
reversed situation. Figure 4 illustrates the lateral and angular displacements of the rotor as 
well as the direction of the fluid film forces and moments. 

For the cylindrical geometry, the film thickness and time derivative are written: 


H = H 0 - (e x +ZB)c os0 - (e y -ZA) sin0 


dH 

dt 


( 



f 



+ Z ™ 

COS0 - 

de y_ z dA 

dt 

dt J 


{ dt 

dt) 


sin0 


( 1 ) 


where H 0 , an arbitrary function of film coordinates, represents the film thickness distribution 
for a rotor that is aligned and centered with the housing. A and B represent the angles of 
rotor rotation about the x and y axes, respectively, while e x and e y represent the components 
of rotor eccentricity at the seal mid-length. The former are referred to as the angular 
displacements and the latter as the lateral displacements. 

For the face geometry, 


H = H 0 + e z - Br cos0 + Arsin0 


dH 

dt 


de 7 dB dA 

— £ - _rcos0+ — rsin0 
dt dt dt 


( 2 ) 


where r is the local radius and e z represents an axial eccentricity of the rotor away from the 
housing. 


One innovative feature of the code is the ability to specify recessed regions in the field with 
bi-directional tapers, as shown on Figure 5, by specifying the depths at the comer points of 
the region and using linear variations between the comer points. 
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901376 


Figure 1. Cylindrical Seal Geometry 
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95TR34-V4 


Figure 3. Face seal geometry. 
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Film Thickness 



Figure 5. Arbitrary Film Thickness Specification 
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2.1 Governing Equations 

The equations governing the flow of incompressible fluids in thin films are obtained [4,5,6] 
by integrating the Navier-Stokes momentum and continuity equations across the film*: 

(f,Re^f b Re b ) n _ _ H z dP ^ (Re^LU Re b f b U b ) 

2 Mr ae 

( f i R ej+f b Re b ) v _ H 2 dP 

2 m h ac 


2 

(3) 


i a 
r ae 


(UH) + 


j__a 

Rd c 


(\/H) + 


dH 

dt 


= 0 


(4) 


where fj and f b are the friction factors relative to the housing and journal surfaces, 
respectively, and are functions of the Reynolds numbers relative to these surfaces as well as 
their roughness. They are given by: 


Re,= £Hij(U-U l Y + V 2 
M 

where i=j,b, and 



1000 ( laminai ) 

(1-35 2 + 25 3 ) +//( 35 2 -2f), 1 000 <Re, <3000 
f‘ , Re j i. 3000 (turbulent) 


(5) 


( 6 ) 


Re r 1000 
2000 


f'= 0.001375 1 


+ 


10 4 e f 

H 


+ 10 6 3 

2 Re, 

1 y 


(7) 


bilm and Ji/m (hickness are used to mean the gap of lubricant separating the rotor and housing. 
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The friction factor for turbulent flow through pipes, f\ in equation (7) uses the curve-fit 
obtained by Nelson [7] to Moody's data. The transition from laminar to turbulent flow is 
obtained using a cubic polynomial which matches values and slopes at both ends, as reflected 
by equation (6). Figure 6 is a plot of the friction factor versus Reynolds number and surface 
roughness, while Figure 7 is an enlargement showing the detail of the transition region. 

Under laminar flow the friction factors are equal to 12/Re and the momentum equations can 
be solved explicitly for the velocities in terms of the pressure gradients: 

U= _}2H^dP + + <* 

\JR 90 2 

( 8 ) 

v= _ 12 H 2 dP 
VR 9{ 

There is no effect of roughness if the flow is laminar. 


Lubrication Background 

In the classical theory of lubrication, when the housing is stationary and the rotor wall velocity is 
Uj=Cor, the fluid velocity components are expressed explicitly in terms of the pressure gradients: 


,, = _H 2 G,0p ^ cor v __H 2 G zdp (9) 

1 2(J r 90 2’ 12pft9{ 


where G x and G z are turbulence coefficients [8] which become unity in the laminar regime. 
Substituting these velocity components into the continuity equation, results in the classical 
Reynold's equation: 


i a 


+ 1 3 

\ h 3 g aP l 

r 2 ae 

x SO J 

R 2 d f 

acj 


= 6pco^ + 12p aH 


ae 


dt 


( 10 ) 


Boundary Conditions 

Boundary conditions on the film pressure distribution consist of prescribing either the pressure 
at the boundaries of the film, the flow normal to these boundaries, or a relation between these 
two quantities. 

At the circumferential ends of the seal surface model, either the pressures are prescribed: 

P= 0 at 0 = 0 S and P= 0 at 0 = 0 e , 
or periodic boundary conditions exists: 

p( e = e s ) = p( e = 0 e ) and u( e = e s ) = u( 0 = e e ) 

Periodic boundary conditions are used, for example, for a 360° seal, where 0 c =O s +27t. 
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Friction factor f 



Figure 6. Friction 
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Friction factor, f 


Fig.1 Transition friction factor 



Figure 7. Detail of Friction Factor in Transition Region 
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For a cylindrical seal the pressure/flow relationship is prescribed by: 

P= P, - K e l / 2 pV n 2 and P= P r - K e l / 2 pV n 2 

at the left and right ends of the seal surface model, respectively. When a symmetry boundary is 
present at the seal mid-length, only the left half (-V 2 L < z < 0) of the seal is modeled and the 
right end relationship is replaced by a zero axial velocity, V(O,0)=O. For a face seal, the same 
relationships are used at the inner and outer radii, although no symmetry boundary is possible. 

P= ^ - K e */2pV n 2 and P= P r - K e ‘/ 2 pV n 2 

At pocket boundaries, the relationship of film pressure to flow is written: 

P= P p - K e ‘/ 2 pV n 2 . 

In all of the above relationships, V n is the flow velocity at the entrance to the film, normal to 
the pressurized boundary and given by: 




V-h, V-n>0 
0, Vh <. 0 


V = Ue z + Ve Q 


(11) 


No pressure drop exists in the case of reverse flow (i.e., flow into the pressurized boundary). 


External Pressurization 

The pressure drop across the orifice supplying the pocket is given by: 


P p = sgn(Q r )^ 


Q, 

A 0 C d 




( 12 ) 


where A 0 is the orifice area, C d is the discharge coefficient and the flow Q r is obtained by 
satisfying continuity over the pocket volume: 


Qr 


§H V h dS + 



(13) 


where A p is the pocket area, S p is its perimeter. The contribution of V*n to this last equation 
may be positive or negative. 
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Dimensionless Variables 


Using the following transformation to dimensionless variables, 


b= B (C 3 /12 jaR 4 ) 
f= F/(P„R 2 ) 
h= H/C 
k= K (C/P 0 R 2 ) 
m= M/(P C R 3 ) 


T= t (C 2 P c /12 gR 2 ) 
A„= 6 gU b R/(C 2 P 0 ) 
Aj= 6 gU,R/(C 2 P 0 ) 
e= e/C 
a= A (R/C) 


P= P/P„ 

q= Q r (12 g/P„C 3 ) 
u= U (12 gR/C 2 P„) 
v= V (12 |oR/C 2 P 0 ) 

^ = Z/R for cylindrical geometry 
= r/R for face geometry 


P= B (R/C) 

Re*= ph 3 Vp/ p 2 
Re c '= pC 3 P,/(R g 2 ) 

A = pC 6 P 0 /(288 A 0 2 C d 2 g 2 ) 
= (Re 0 */288) (C 3 R/A 0 2 C d 2 ) 
A = K e (Re 0 *C/288 R), 


the dimensionless film thickness is now written for the cylindrical geometry: 

h = h a - (e x +zp)cos0 - (e y -za)sin0 

. ...1 

sin0 



f 

A 


C \ 

3«y_,0« 

dh _ _ 

3e x 

+ zl£ 

COS0 - 

dx 


0T, 


, dx dx J 


and for the face geometry: 


h = h Q + e z - p C cos0 + a C sin0 

dh de z gn g a . 

— = — - - _li£cos0+ — fsin0 

dx ox ox dx 


Equations (3), (4) and (5) become: 


(f^e^Re,) u __ _ 12/j2 |g + 


(f,Re i+ f b Re b ) 


00 

v = -1 2h 2 — 
dz 
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-|(uh) ♦ Mvti) ♦ - 0 

00 dz dx 


( 17 ) 


Re' h I : 

Re r ^-^(u-2 A,.) 2 + v 2 , i=j,b 


( 18 ) 


Equations (6) and (7) remained unaltered, as they were already dimensionless. 

The dimensionless form of the boundary conditions now become, at the circumferential 
ends, either 

P(t e.) = 0 and p(£, e e ) = 0 
or 

P(C» 0.) = P(C» 0e) and U (C 0 S ) = U(^, 0 e ) 

when periodic boundary conditions are present. At the left end, the dimensionless boundary 
conditions become 

p(-L/D,0) = Pl - A e v n 2 

and at the right end, either 

p(L/D,0) = p r - A e v n 2 
or 

v(O,0) = 0. 

At pocket boundaries 

P(z,6) = p p - A e v n 2 

where: 


^ < 


Vh, 

0, 


v = ue. 


V-n> 0 

+ ^e 


Equations (12) and (13) governing the external pressurization become: 

Ps - Pp = sgn(q r )A r q? 


( 19 ) 


( 20 ) 
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q r = 


fh 

s „ 



— dddz 
dx 


( 21 ) 


2.2 Solution of Film Pressures 

Discretization of the seal surface is done by using a rectangular grid, with M lines in the 
meridional direction and N lines in the circumferential direction. The grid lines are separated 
by variable increments. The pressure distribution is represented by discrete values at the grid 
points located at the intersections of the grid lines. There must be grid lines coincident with 
the boundaries of the seal surface and with the pocket boundaries. Using the cell method [9], 
a control area or cell is centered at each grid point and extending half way to the neighboring 
grid lines, as shown by the shaded area in Figure 8. The grid points are noted by the solid 
circles and have grid coordinates i,j. The film thickness is evaluated at the comers of the cell 
(denoted by the shaded circles marked h t , h 2 , h 3 , and h 4 ) located at the geometric centers of 
the rectangles formed by the grid lines. This staggered configuration allows a discontinuous 
film thickness to be treated, as occurs, for example in a seal with a Rayleigh-step. 
Circumferential and meridional components of velocity are also associated with each of the 
four cell comers. 

Using the divergence theorem, the continuity equation may be integrated over the cell to give: 

-j hvn dS = j dA (22) 

S c K 

where A c and S c are the cell area and perimeters, respectively. The left hand side of the 
above equation is the sum of the flows out of the cell while the right hand side is the rate of 
change of the cell volume. The finite-difference form of this equation is: 


F u = ) + -u 3 h 3 ) + 

+ ~ v * h *) + ~v 3 h 3 ) - (23) 

" J yjf ( Az / + Az mX A0 ( + A 0 m) = 0 


where F g is the error in satisfying continuity of flow in the cell centered at i,j. Although the 
time rate of change of film thickness has been evaluated at the center of the cell, it could 
have alternatively been evaluated at each of the four cell comers. 
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When the grid point falls on a pressurized boundary, such as a pocket or seal end, the film 
pressure error is: 


F u = P b - Pu ~ K max(0 , v n f = 0 


v n = 


_ S, 


S b h;,i 


(24) 


where p b is the dimensionless boundary pressure*, v n is the mean velocity of the flow that 
crosses the portion of the boundary perimeter that intersects the cell, and Z (j represents the 
sum of the appropriate terms in equation (23) contributing to the cell flow. Figure 9 shows 
an example of the cell i,j located at the right bottom comer of a pocket. In this case, the 
mean velocity would be evaluated as: 




A z, 


L (uA) 


A z 




AG; 


A 0 1 

_ 7-1 


2 ( u 2 h z -“A) + -v 2 h 2 ) - ~ 


0/7, v (Az,+Az H )A 0 / + Az m A6 h 


(A 6,-1 + A Z/) hjj 

0T 4 


2 


(25) 


Equations (23) and (24) represent the finite-difference form of the continuity equation that 
must be solved for the pressures. The eight components of velocity used in these equations 
are functions of the nine pressures at or neighboring grid point i,j, and are evaluated as 
described in section 2.3. Following the procedure described in Reference 10, these highly 
nonlinear equations can be solved using the Newton-Raphson iteration method [11]. The 
procedure is started with an initially guessed or previously calculated pressure distribution, p ;J . 
The error function Fy is then linearized about this guess in order to obtain a better 
approximation to the pressures p kj new : 





(26) 


where a forward difference or a central difference may optionally be used to numerically 
evaluate the partial derivatives. Pressures without the superscript new relate to the previous 
or "old" approximation. If we introduce the column vector {Pj new } as the M new pressures at 
the jth column of grid points. Equation (26) may be written: 

[C i ]{Pj r>ew } + [E'HPiT} + [D'KPj." 6 ”} ={R'} , < 27 ) 


*P/P 0 , P/P 0 or P/P 0 . 
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95TR34-V4 

Figure 9. Example of Cell at Corner of Pocket 
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where [C j ], [E j ] and [D j ] are tri-diagonal matrices whose interior elements are: 


C 1 

v -'i,i+k ~ 


0F 1; 


ap 


E.i = 


> *— i,i+k 


i+k,j 


9P: 


D 1 

> '-'i,i+k 


aF ;; 


i+k,j-1 


ap. 


, k = -1,0, 1 ; i = 2 M -1 


i +k,j +1 


The interior elements of the column vector {R j } are: 


- (^U+kPi+k.j + J+kPi+k,j-1 + J+k Pi-tk.j+1 ) - ^"ij 

k = -1 


The set of linear equations (27) that result for the next guess of pressure distribution is in a 
form suitable for solution by the column method which is described in detail in References 9 
and 12. This method makes use of the banded nature of the equations in order to minimize 
computer time. 


2.3 Solution of Flow Velocity 

The momentum equations (16) are used in order to evaluate the velocity components from the 
pressure gradients. These equations may be rewritten in the generic form: 

e /A +fl VA) = o. 

(28) 


where the Reynolds numbers used to evaluate the friction factors are based on the magnitude 
of the local fluid velocity relative to each surface: 


dp 

00 

0p 

0Z 


fU,V\ 


,u t v\ 


( R 


f,Re,+f h Re h „ a n 

= — — !—Z — i u + 12/) 2 — 

2 ae 

f : Re i +f h Re tl „ 

= — — — i v + 12 h 2 ^R. = 0, 
2 dz 


Re r - 2^/-(7-2Ay + , 2 , 

Re*h I 

R e b -~i (u~2A b f + , 


(29) 


The dependence of the friction factors on velocity components orthogonal to each momentum 
direction couples the two momentum equations. Figure 10 is a schematic of the rectangular 
region between meridional grid lines i and i+1 and circumferential grid lines j and j+1. In 
order to preserve continuity, it is essential that the same equation be used to evaluate the 
velocity components for adjacent cells. That is, the velocity u, out of the shaded cell centered 
at ij must have the same value as the velocity u 4 into the cell centered at i,j+l (see Figure 8). 
This value is designated as u' in the figure. Similarly, the velocity v, out of the cell ij must 
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Figure 10. Schematic of Rectangular Region Between Grid Lines 
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be the same as v 2 into the cell at i+l,j, and is designated as v\ This is achieved by using the 
average of the two corresponding orthogonal components. Thus, the component u' is 
determined by the u-momentum equation: 



Pij+rPij 




v+v + 
2 


= 0 


while the component v is determined by the v-momentum equation: 



Pmj-P., 
A z. 


LT+U 


V 


= 0 


Similarly, u + and v + are determined by: 


( 30 ) 


( 31 ) 


Pi+\ ,;'+1 P MJ + V +V + 
> ! ' 



( 32 ) 


Equations (30), (31) and (32) are four coupled equations that determine the velocity 
components from the four pressures at the corners of the rectangle between grid lines and 
must be solved simultaneously. This is done using an inner Newton-Raphson iteration loop. 
By performing the differentiation of the error functions (G u , G v , ...) with respect to the four 
unknown velocities, analytically instead of numerically, significant computer time is saved. If 
the velocities have not been previously calculated initial guesses may be obtained from equations 
(8) assuming laminar flow. Once the iterations for the velocities have converged, their values are 
saved to provide a good starting guess for the next time they must be calculated. 

One simplification is possible by assuming that the friction factors are constant within the 
rectangular region and the Reynolds numbers are based on the averaged flow velocity 
components, Vi(u'+u + ) and V 2 (v +v + ). Although this does not uncouple the four equations, it 
requires less number of evaluations of the square root in equation (18). Since this assumption 
saves some computer time without introducing significant errors, it was chosen as the default 
program option. However, occasionally when the grid is not very fine and the pressure 
gradients vary rapidly, the iterations will diverge and the more rigorous formulation, which 
uses distinct friction factors for each of the four momentum equations, should be used. 

If the surfaces are smooth and the housing is stationary so that the continuity equation takes 
the form of equation (10), the simpler formulation described in detail in Reference 10 may be 
used, resulting in significant reduction in computer time. 
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2.4 Fluid Film Load, Moment and Torque 

The forces and moments on the rotor generated by the fluid film pressure distribution are 
obtained by integration of the pressure distribution over the cylindrical seal surface: 





F x 


COS0 

fy 

l e 9 

sin0 


'-u H 
~ L 0 S 

-Zsin0 

M x 

kl 


Zcos0 


(RdQ dZ 


The dimensionless form of this equation is written: 


f. 

L 

' ' 

COS 0 

. f y 

D 8 S 

" C C 

sin 0 


= \ \ P' 



_ L 0 S 

-zsin 0 

m y 

D 

z cos 0 


\dQ dz 


( 33 ) 


(34) 


For the face geometry, the differential of forces and moments are given by 



[Fj 


' > 

1 

d < 


► = p' 

-rsin0 


KJ 


rcos0 


(35) 


where the pressure p is a function of r and 0 and 

dA = r dr dd (36) 

The pressure in the rectangular region 0j<0<0 j+1 , R i <r<R i+1 can be expressed most accurately 
as a bilinear function of the coordinates as: 

p(r,d) = a 1 + a 2 r + a 3 0 + a 4 rd (37) 
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where the coefficients are functions of the values of the pressure at the four comers of the 
rectangle: 


a _ (P/.,A. -P»u R ,)8j. , 

1 (« M - R ) (0;.. - 0, ) 
a _ Pij + 1 R+1 Pi+^,l+^ R i Pij R i*\ + Pi*1,j R i 

2 (A w - R) ( 9 /*i - 0/ ) 

a _ (fi'+lj ) 0 .- 1 {Pi* 1./+1 Pi, j* 1 j 9 ; 

3 " («M - *,) (9/*i - 0; ) 

q _ Ph1j*\ Pi,'}*} A+1,/ + Pij 

4 " ^ (9,, - 9, ) 


( 38 ) 


The contribution to the fluid film forces and moments corresponding to the rectangle is then 
obtained by integrating equation (37): 


to obtain: 



> 

F, 

®/«-i Fi M 

a t r + a 2 r 2 + a 3 rtJ + a 4 r 2 0 

A’ 

M x 

= J J • 

(-a^ 2 - a 2 r 3 - a 3 r 2 0 - a 4 r 3 0)sin0> 


KJ 

9; R, 

(a 1 r 2 + a 2 r 3 + a 3 r 2 0 + a 4 r 3 0)cos0 


drdd 


, ft 


'f; 


’ SrS? 


SUe 

SfflSa 

-M x 

► ~ 

_cx 2 cx° 


_qj 2 cx 1 

■ofl *o s 

-cx 3 cx 1 
•Or •Os 

M y. 


ex 2 cx 3 


cx 2 cx 1 
•Or «o c 

cx 3 cx 1 
•Or 


P4 


( 39 ) 


( 40 ) 
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where 


a" = f 0 n cos0 c/0 = [ cos n 0 - 0 n sin0 ] , 

*> 0, 

0 , 

^ 0 ^ 

Ss = J 0 n sin0 c/0 = [ sin n 0 - 0 n cos0 ] , 
0 / 9; 


a 


n _ 

0 = 


J 0,1 de 


,n+1 


0 " 1 o' 
7+1 Oy 

nTT 


n+1 


Sp = J r n dr = r n dr = 


n n+ 1 q' 
rij 

n + 1 


n+1 


For the cylindrical geometry, 


' 



F x 


-COS0 

■ F y 

M x 

■4 

' =P 

-sin0 

zsin0 

KJ 


-ZCOS0 


RdQDZ 


( 41 ) 


(42) 


The pressure in the rectangle (0j<0<0 j+1 ,z i <z<z i+1 ) can be expressed as a linear function of the 
coordinates : 


p(z,0) = a 1 + a 2 z + a 3 0 + a 4 z0 


(43) 
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where the coefficients are functions of the values of the pressure at the four comers of the 
rectangle: 


a _ Pf+1 ,/ (PjJ+\ Z i + 1 z /) 0 / 

( z m - 4 ( 0 m - 0 , ) 

a = Pi ^ z / Pjj z /+i + Pj+ij^j 

( Z /+1 ” Z /) (®/+1 ” ) 

a _ Pi,j ) 0 y +1 (P/+i,/+i ~ P/,/ 4 i )Qy 

( z m " ( 0 m - 0 , ) 

a _ Pmj+i Pjj*\ Pmj + Pj,j 

Ki - 4 ( 0 /*< - 0 y ) 


Substituting equations (43) & (44) allows integration of (42) over the rectangle: 
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(45) 


(46) 


( 47 ) 
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The differential of torque transmitted from the housing to the rotor is given by the cross 
product of the position vector r and the shear traction vector acting on the housing T: 


T = Te z = ffp x t d A 

A, 

= R JJe f x ( T *^r) dA 


( 48 ) 


T = 


P„R‘ 


ff 9 p _ f M u-2 Aj) - f b R„(U-2A b ) 
JJ 1 aa 72/7 1 


2 c 0 ae 

O A, 


For laminar regime, f J Re J = f b Re b = 12, and the equation simplifies to: 


T = 


p ofl 2 rc, h dp _ 


2 C 


^ Mdz 


o A, 


ae 3/7 


(49) 


The power loss due to the difference in velocities across the two surfaces is obtained by the 
dot product of this torque with the relative velocity: 


P o* 2 ffru9p 
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(50) 
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2.5 Stiffness and Damping Coefficients 

Defining W to be a generalized vector of forces and moments generated by the fluid film 
pressure and r to be a generalized vector of lateral and angular displacements: 


f, 

A 
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e„ 

► - a 

e y 
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y 

r= — 
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yj 


a 

iPj 

dx 
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the matrices of stiffness and damping coefficients can be written: 


( 51 ) 


1 A 


b . - 

’ df , i 


(52) 


where the subscripts i and j range over x, y, a and [3. These coefficients are evaluated by 
numerical differentiation of W, using a forward difference. For example: 



F y (e x , €y, qtg, p) - F e y , cc, p) 

5 


(53) 


2.6 Solution of Rotor Position and Pocket Pressures 

If the rotor position is specified, equation (51) is used to solve for the fluid film forces and 
moments in terms of the calculated pressure field. Similarly, if the pocket pressures are 
specified, equation (12) is used to solve for the orifice size in terms of the supply pressure 
and calculated pocket flow. 

On the other hand, if externally applied loads and moments on the rotor (f , f , m vor and m vo 
are specified they must be balanced by the fluid film forces to maintain static equilibrium. 
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Similarly, once the orifice size is specified, equation (12) must be satisfied by the pressure in 
each pocket. The global set of equations that must be satisfied by the rotor displacements and 
pocket pressures are: 


m - a, 

w - - f y g 

m M) = - m x 9 
m y (r) = -m yg 


Ps - Ppi = sng(q r1 ) A r , (q r ,f, for pocket 1 , 

I 

Ps - Pp2 = sn 9(qJ A^(q r2 ) 2 ’ P° cket 2, eto. 


(54) 


The vector r can now be redefined to include the pocket pressures and a generalized vector 
of errors in forces, moments and pocket pressures W e can be defined: 
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(55) 


Solution of the global equations is performed by Newton-Raphson iterations, as follows: 


W + 


dW e 

dr 


( r nm - r ) = 0 


(56) 


where, as before, the superscript new indicates the newer values of the vector r. 
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3.0 SAMPLE PROBLEMS FOR CODE ICYL 


A number of sample problems were completed to demonstrate the behavior and various features 
of the computer code. They are intended primarily for illustration and do not necessarily 
represent recommended seal designs. Table 1 summarizes the mesh size and approximate 
execution times (on a 20 MHz 386 PC), as well as the specified and computed variables. 

The complete input file is included at the top of the output file for each case. The filenames 
used in the samples all used the prefix icyl. For example, EX1 used /name = icylEXl. 

Samples EX1, EX2, and EX3 were selected with a coarse (5 x 11) mesh covering a 90° sector 
in order to demonstrate the use of pressurized pockets and iterations for rotor position within 
a reasonable execution time. A pocket with a supply pressure of 100 psi was centered on the 
seal sector modeled. 

Sample EX1 contained two cases. In the first case, the pocket pressure was specified as 50 psi, 
resulting in an orifice diameter of 0.0137 in. calculated by the program. Both components of the 
resulting fluid film force are equal and the moments are zero, as would be expected at the 
concentric position. In the second case, the rotor was moved with eccentricity ratio of £ x = 0.1 
and given a misalignment ratio of |3 = 0.1 about the y-axis. With the value of orifice 
diameter already assigned from the first case, the pocket pressure and forces rise slightly, 
generating non-zero moments. 

In sample EX2, external forces and moments equal to the negative of those resulting in EX1 
were specified, in order to have the program iterate for the rotor radial and angular positions. 
Five unknowns, the four displacement components (e x , £ y , a, and P) as well as the pocket 
pressure, are iterated for simultaneously. The initial pressure distribution was read from the 
converged values of EX1. Since the iteration was begun at the concentric position where the 
orifice was sized, the pocket flow error was zero and increased when the rotor was moved in 
the first iteration, causing the ran to abort. When the limit on diverging iterations was 
increased to 2, the iterations converged in only 3 iterations to within a small error of values 
expected (£ x = 0.1, P = 0.1). 

In sample EX3, an axial taper of ±30% of the clearance from end to end was superimposed. 

In this sample, the program was asked to find the angular rotor position such that no external 
moments were required, while the rotor eccentricity was varied, using £ x = 0.1, 0.3 and 0.5 
for the first, second, and third cases, respectively. The results show that as the eccentricity is 
increased in the x-direction, the rotor twists about the y-axis in order for the moments 
generated by the film to be zero, resulting in p = 0.050 and 0.17 at £ x = 0.1 and 0.3, 
respectively. The case of £ x = 0.5 resulted in a negative film thickness with the appropriate 
error message and recommended user action: 

• Reduce the specified applied forces/moments 

• Reduce the specified eccentricity/misalignment. 

The resulting film thickness and pressure distributions for sample EX3 are shown in Figures 
11 and 12, respectively. 


NAS A/CR— 2004-2 1 3 1 99/VOL4 


31 



Table 1. Summary of Sample Cases for Code ICYL 


Case 

Mesh 
Si ze 

I SYM 

Variable 

Found Specified 

App.exec 

time" 

Features 

EX1 

5x11 

0 

dor ’'f.Pp<x* 

Pp<x**^x»dorif 

4.6 min 

1-pocket 

EX2 

5x11 

0 

^xt P»Ppock Pyf ^ 

11.6 min 

1-pocket 

EX3 

5x11 

0 


M xf N y ,dorif 

6.2 min 

Tapered pocket 

F3 

9x61 

1 ' 

- 

all 

29 min 

Raleigh-step 

F4 

7x61 

0 

- 

all 

7.8 min 

Axial taper 

F9 

5x73 

1 

K,8 

all (3 preloads) 

1.6 hrs 

3- lobe 

11 

5x61 

1 

dor if 

P(»<* 

7 min 

4-pocket 

12 

5x61 

1 

P if 

r podtf ^ 

c x ,dorif 

1 .8 hrs 

4-pocket 

13 

5x61 

1 

^xt P pock 

F xl F y ,dorif 

1.9 hrs 

4 -pocket 

14 

9x61 

0 

KfB, P,*** 

« x ,a,dorif 

7.7 hrs 

4-pocket 

15 

11x61 

0 

dor i f 

Ppoc*. 

5.2 min 

8-pocket 

16 

11x61 

0 

Ppw* 

dorif ,f x ,a 

3.1 hrs 

8-pocket 

015 

5x31 

0 

K,B 

all 

1 hr 45 sec Roughness 


K,B indicate whether stiffness, damping coefficients were requested, 
on an IBM PS/2 model 70 (386 20-Mhz) computer. 
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ICYL: ICYLEX3 FILM TH I CHNESS DISTRIBUTION 

90 DEG SECTOR: CALCULATE ANGLULAR POSITION FOR ZERO MO 
RAD I US= 1,0000 SPEED- 1000.0 

LENGTH= 4.0000 

CLEARANCE^ 1000QE - 02 
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Figure 11. Film Thickness Distribution for Sample EX 3 
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Sample F3 shows a 120° sector with a Raliegh step of linearly varying depth. The resulting 
film thickness and pressure distributions are shown in Figures 13 and 14, respectively. 

Sample F4 shows a 120° sector with an axial taper in the right half (4<i<7) of 0.001 in. 

Since two fewer intervals were used in the axial direction than in the previous cases, and 
since half as many iterations were required for the pressure distributions, the execution time 
was reduced from about 29 to 8 minutes. The resulting film thickness and pressure 
distributions are shown in Figures 15 and 16, respectively. 

Sample F9 is that of a full 360° seal with three 60° lobes. The dynamic coefficients were 
requested as the preload was increased from 0.1 to 0.3 in the middle case, to 0.8 at the last 
case. The resulting film thickness and pressure distributions are shown in Figures 17 and 18, 
respectively. 

Samples II through 114 are more realistic models since they represent the full circumference. 
In sample II, the program takes about 7 minutes to calculate the orifice diameter for given 
pocket pressures. The resulting pressure distribution is shown in Figure 19. Sample 12 is the 
same as II except that the eccentricity and orifice diameter were prescribed, requiring the 
program to solve for the four pocket pressures. The resulting film thickness and pressure 
distributions are shown in Figures 20 and 21, respectively. In sample 13 , the orifice diameter 
as well as the radial force were prescribed, requiring the program to solve for the radial 
position as well as the pocket pressures. 

Sample 14 shows the dynamic increase in execution time with the number of axial grid lines, 
M. Sample II is a model of only half of the seal, symmetry imposed, at the concentric 
position with the pocket pressures specified. This run executes in less than 7 minutes in spite 
of the 5x61 mesh. Sample 14 is a model of the full axial length (IS YM = 0), with an 1 1 x 
61 mesh, in which non-zero £ x , a and orifice size are specified and all 32 dynamic coefficients 
are requested. This run took 7.7 hours to execute. Calculations for each of these coefficients 
require convergence of the outer iteration loop with four unknown pocket pressures. 

Samples 15 and 16 are models representing the full circumference and length with two rows of 4 
pockets. The orifice size is calculated in the concentric aligned position in 15 while 16 calculates 
what happens when the rotor is displaced to the £ x = 0.4 position and rotated about the x-axis by 
a = 0.4. For 15, the resulting pressure distribution is shown in Figure 22. For 16, the film 
thickness and pressure distributions are shown in Figure 23 and Figure 24, respectively. 

Sample 015 is a case of a plain cylindrical seal with increasing housing wall roughness. The 
wall roughness (ROUGHB) was varied from 1 x 10 6 to 1 x 10 3 in. in a logarithmic scale. 
This input was used to generate the top curve of critical mass versus roughness shown in 
Figure 25. The stabilizing effect of housing roughness is more pronounced at the higher 
pressures due to the increased effect of inlet inertia. 
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ICYL: ! CYL F3 FILM THICHNESS DISTRIBUTION 

RALEIGH STEP WITH VARIABLE DEPTH 

RADIUS^ 4.0000 SPEED= 1000.0 

LENGTH^ 2.0000 

CLEARANCE^ . 10000E-02 
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Figure 13. Film Thickness Distribution for Sample F3 
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ICYL: ICYLF3 PRESSURE DISTRIBUTION 

RALEIGH STEP WITH VARIABLE DEPTH 

RADIUS^ 4.0000 SPEEDS 1 000. Q 

L ENGTH= 2.0000 

CLEARANCE^ . 1 00 DOE - 02 
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Figure 14. Pressure Distribution for Sample F3 
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ICYL: ICYLH FILM THICHNESS DISTRIBUTION 

TEST OF AXIAL FILM TAPER 

RAD I US= 4.0000 SPEEDS 1 000.0 

L ENGTH= 2.0000 

CLEARANCE= . 1 000DE - 02 



95TR34-V4 


Figure 15. Film Thickness Distribution for Sample F4 
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Figure 16. Pressure Distribution for Sample F4 
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DELTA(1 ,7) 
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ICYL ICYLF9 

PRESSURE DISTR 

1 BUT ION 


3-LOBE SEAL 

PRELOAD^ 0.8 



RADIUS^ 

4 . 0 0 D 0 

SPEEDS 

1000.0 

LENGTH= 

2.0000 



CLEARANCE= 

. 1 0OOOE - 02 
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Figure 18. Pressure Distribution for Sample F9 
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Figure 19. Pressure Distribution for Sample II 
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Figure 20. Film Thickness Distribution for Sample 12 
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Figure 21. Pressure Distribution for Sample 12 
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ICYL: icy! i 5 


PRESSURE DISTRIBUTION 


TEST OF ITERAT 

IONS FOR 

POCKET PRESSURES 


RADIUS^ 1 

0000 

SPEED= 

3600.0 

LENGTH= 1 

0000 



CLEARANCE= 

. 1 0OOOE 

-02 
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Figure 22. Pressure Distribution for Sample 15 
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Figure 23. Film Thickness Distribution for Sample 16 


NAS A/CR— 2004-2 1 3 1 99/VOL4 


46 



ICYL: ICYLI6 PRESSURE DISTRIBUTION 
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Figure 24. Pressure Distribution for Sample 16 
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4.0 SAMPLE PROBLEMS FOR CODE IFACE 


A number of sample problems were completed to demonstrate the behavior and various features 
of the computer code. They are intended primarily for illustration and do not necessarily 
represent recommended seal designs. Table 2 summarizes the mesh size, approximate execution 
times (on a 33 MHz 486 PC) as well as the specified and computed variables.* 

Sample 1 uses a coarse mesh, variable in both directions, to represent a 50° sector with a film 
thickness that is tapered circumferentially by 1 mil. This sample demonstrates how different 
input variables control the film thickness distribution. The left bottom comer point increases 
the film from the reference clearance to 2 mils. DELTA(3,1) (see Figure 5) then decreases 
the film thickness by 1 mil over the circumferential extent. The axial eccentricity ratio (EX= 
-0.1) decreases the clearance distribution by 10% of CREF. The model is located in the first 
quadrant, so the resulting pressure distribution produces a positive moment about the x axis 
and a negative one about the y axis. The resulting film thickness distribution is shown in 
Figure 26 and the pressure is shown in Figure 27. 

Sample 2 A specifies a geometry that is repeated circumferentially every 120°. A zero 
pressure is specified at the first two grid lines of each pad while the rest of the pad's film 
thickness is tapered. The prescribed external force and moments were obtained from the 
negative of the results of a similar case at e z =0.2, a=0.5 and (3=0. In sample 2A, the 
program iterates on the position, starting from zero guess values, since Ex, ALFA or BETA 
were not prescribed, and converges on e z = 0.2, a = 0.500002, and a near zero (3 in 4 outer 
loop iterations. The resulting film thickness distribution is shown in Figure 28 and the 
pressure is shown in Figure 29. 

Samples 3A, 3B and 3C are used to illustrate the steps that would be used in calculating the 
performance of a hydrostatically supported, 4-pocket face seal. A variable mesh is used to 
describe the geometry of l A of the circumference between 0 = 0 and 90° from the x-axis. One 
pocket with a supply pressure of 1000 psi is centered on the section modeled. Pressures of 100 
and 200 psi are prescribed at the inner and outer radii. Periodic circumferential boundary 
conditions and 4 pads are specified to generate the geometry in the remaining three quadrants. 

In sample 3 A, sizing of the orifice diameter is performed with pocket pressures at 500 psi and 
an aligned rotor. Although the l A model corresponding to one pocket would be sufficient for 
the orifice calculation, the full model was required in this case in order to calculate the 
angular mode stiffness and damping coefficients. The output indicates an orifice diameter of 
0.069 in., a total flowrate of 9 in. 3 /s, and a fluid film axial force of 7,376 lbf. The direct axial 
stiffness is K zz =12.7xl0 6 lb/in. and the direct angular stiffness and damping coefficients are 
33.2xl0 6 lb-in./rad and 29.6xl0 3 lb-in./rad. 


K and B indicate calculation of stiffness and damping coefficients. 
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Table 2 . Summary of Sample Cases 


Case 

Mesh size 
MxN 

Variables 

specified 

Variables 

calculated 

NPADS 

run time 
(sec) 

features 

1 

5x11 

EX, ALFA 



3.75 

variable grid, DELTA(l.l) 

2A 

5x31 

FXG, MXG, 
MYG 

EX, ALFA, 
BETA 

3 

125 

prescribed force & moments, 
DELTA(5,1) 

3A 

7x41 

PPOCK 

DORIF, K,B 

4 

404 

4-pocket, calculation of orifice & 
coefficients 

3B 

7x41 

DORIF, 
EX, BETA 

PPOCK 

4 

352 

4-pocket, prescribed 
displacements 

3C 

7x41 

DORIF, FXG, 
MXG, MYG 

PPOCK, EX, 
ALFA,BETA 

4 

566 

4-pocket, prescribed force & 
moments, pressures read 

13 

9x37 


K, B 

4 

1093 

preloaded pads, roughness 
multiple cases, DELTA(4,1) 

17 

9x65 

ALFA 


8 

103 

8 Rayleigh steps 

15A 

10x61 

PPOCK 

DORIF, K,B 

4 

2115 

4-pocket with XKE=1 

I5B 

10x61 

DORIF, EX, 
MXG, MYG 

ALFA,BETA, 

PPOCK 

4 

1368 

4-pocket finding angular position, 
pressures read 


95TR34-V4 


NAS A/CR— 2004-2 1 3 1 99/VOL4 


50 



























































I f A C r : 1 1 A C I I f 1 1 U I II I C II H li 5 S D I S I R I II U I I 0 N 

I F A Cf 1 I* l AIN IIDGI I! Ill VARIAIII! GRID 

RADIUS* 1. 0 0 D D S I' I L I) 1 0 6 0.0 

KNGIII 3. 0 0116 

C l f A R A N C l - 0 . 1 0 0 0 0 1 0 ? 



95TR34-V4 


Figure 26. Film Thickness Distribution for Sample 1 
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Figure 27. Pressure Distribution for Sample l 
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Figure 28. Film Thickness Distribution for Sample 2 A 
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Figure 29. Pressure Distribution for Sample 2 A 
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In sample 3B, the fluid film forces are calculated for prescribed orifice diameters of 0.068 in., an 
axial eccentricity ratio* of 20% and a misalignment ratio of 50% about the y-axis (DOREF, EX, 
BETA). The pocket pressures converged in 4 iterations. The pressure increases in pockets 1 and 
4 and decreases in 2 and 3. The fluid film axial force increases to 7,932 lbf, with resulting 
moments of -287 and -5,216 in. -lb about the x and y axes, respectively. The resulting film 
thickness and pressure distributions are shown in Figures 30 and 31, respectively. 

In sample 3C, the external force and moment are prescribed equal to the negative of the fluid 
film forces calculated in sample 3B. The initial guess for pocket and film pressure 
distributions were read from the values saved previously in sample case 3A. Although 
convergence is also achieved in 4 outer loop iterations, it takes more time to execute (572 
versus 360 seconds) because the program must now simultaneously iterate on three position 
variables as well as the four pocket pressures. 

Sample 13 illustrates multiple inputs in the same run. The model consists of 4 preloaded 
pads, each with a preload ratio of 50% and a circumferential offset of 10°. The first input 
prescribes a roughness of 0.02 mils to both surfaces while the second and third inputs have 
smooth housing and rough rotor, and vice-versa. The last input restores smooth surfaces for 
both housing and rotor. The resulting film thickness is shown in Figure 32; the pressure 
distribution is shown in Figure 33. Table 3 illustrates the effects of surface roughness on the 
torque and direct stiffness coefficients. It is noted that roughness of the housing alone can 
increase the stiffnesses by 43 to 52% with less than 18% increase in torque, while roughness 
of the rotor alone decreases the stiffnesses with equivalent increase in torque. 

Sample 17 illustrates a hydroynamically supported face seal. The sealed pressure at the inner 
radius is 1000 psi higher than the outer radius. Eight Rayleigh steps are equally spaced along 
the circumference and are fed from the inner radius by a deep radial groove. A misalignment 
ratio of 50% is specified about the x axis. The resulting film thickness is shown in Figure 34 
and the pressure distribution is shown in Figure 35. 

Samples 15 A and 15B are similar to 3 A and 3C, respectively, but with a finer mesh, which 
takes significantly longer execution times (2116 and 1702 sec). The resulting film thickness 
and pressure distributions are shown in Figure 36 and Figure 37. 


Prescribing EX=0.2 with CREF=0.00l is analogous to prescribing CREF=0.0012 with EX=0.0. 
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Figure 30. Film Thickness Distribution for Sample 3B 
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Figure 31. Pressure Distribution for Sample 3B 
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Figure 32. Film Thickness Distribution for Sample 13 
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Figure 33. Pressure Distribution for Sample 13 
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Table 3. Effect of Roughness on Torque and Direct Stiffnesses 


roughness (mils) 

torque 


K aa - Kpp 

rotor 

housing 

(in.-lb) 

(10 8 lb/in.) 

(10 8 lb-inyrad) 

0.02 

0.02 

2,900 

245 

1,956 

0.02 

0.00 

2,354 

158 

1,153 

0.00 

0.02 

2,366 

263 

2,035 

0.00 

0.00 

2,006 

184 

1,337 
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Figure 34. Film Thickness Distribution for Sample 17 
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Figure 35. Pressure Distribution for Sample 17 


NAS A/CR— 2004-2 1 3 1 99/VOL4 


57 





I i-ALfc : i f ac e 1 j 

RRES5L 

RE DISTRIBUTION 

ir.\CE15S -• ■ POCI 

:et face seal ppe 

SCRIBED EX and MOMENTS 

RADIUS- 4 

0000 

SPEED: 30000. 

LENGTH: 3 

0000 


CLEARANCE = 

0. 30000E-02 




95TR34-V4 


Figure 37. Pressure Distribution for Sample 15B 
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5.0 VERIFICATION 


ICYL has been compared with the results of two other MTI computer codes as well as with 
currently published data. The first comparison was against a generic bearing program with 
many similar capabilities (GBEAR) based on the turbulent lubrication theory of Ng and Pan. 
A second comparison against a laminar bearing program (GASBEAR) was used to verify the 
calculations of moments and angular coefficients. Finally, comparisons were made against 
calculations published by San Andres in Reference 14. 


5.1 Comparison Against Other MTI Codes 

The first of the MTI computer codes is GBEAR which is fully described in Reference 10. 

This program is based on the turbulent lubrication theory of Ng and Pan [13] and does not 
include surface roughness, housing rotation or calculation of misalignment coefficients. It 
includes inertia pressure drop at exit from pockets but not from seal ends. 

Calculations were made with a 90° seal sector at an eccentricity ratio of 0.5 and with a pocket at 
its center with a prescribed pressure ratio of 0.5. Table 4 shows a comparison of pocket flow, 
orifice size, force, and stiffness and damping coefficients. As expected, comparisons of GBEAR 
against ICYL with the same friction model yielded nearly identical results. With the new friction 
model that includes surface roughness effects, ICYL calculates lower torque (-32%), lower pocket 
flow (-13%) and orifice size (-7%), and force components (-6%). Very good agreement in the 
stiffness coefficients (-4%), and slightly higher damping coefficients (+13%) are obtained. 

Other comparisons against GBEAR in the laminar regime and without pockets yielded 
identical results. 

A second MTI computer code with the fluid compressibility turned off (GASBEAR) was used 
to verify the calculation of the 24 stiffness and damping coefficients which involve rotor 
misalignment. GASBEAR was written for use in conjunction with plane journal bearings and 
cylindrical seals and does not treat turbulence or pressurized pockets. The comparison, in the 
laminar regime and with the same finite difference mesh, yielded identical coefficients. 


5.2 Comparison Against Published Data 

A detailed comparison was made of the five-pad hydrostatic bearing discussed by San Andres 
in Reference 13. This high speed hybrid journal bearing operates at relatively high levels of 
pressurization and relatively low viscosity lubricants, in which the effects of pressure-induced 
turbulence become important. Fluid inertia may also be important. Figure 38 is a plot of the 
pressure distribution at the concentric position, while Figure 39 and Figure 40 plot it for 40% 
eccentricity ratio of the journal between pockets and over a pocket center, respectively. 
Reproductions of the corresponding pressure distributions published by San Andres are 
included in the figures for comparison. It is noticed that the size of the pressure drops at the 
pocket exits (i.e., entrance to the film) as well as the general pressure distribution are 
comparable for both analyses. 
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Table 4. Comparison Against GBEAR 



GBEAR 

ICYL 

IFRIC=0 

ICYL 

IFRIC=3 

ICYL 

IFRIC=4 

iJimiii.iiiii.uiM 

25.75 

25.21 

20.931 

22.316 

KRHmmiMi 

0.0833 

0.0820 

0.0752 

0.0776 

Torque (lb-in) 

14.38 

14.32 

8.791 

9.771 

Power (Lb-in/s) 

45,171 

44,971 

27,617 

30,696 

Fx 



3,694 

3,358 

3,352 

3 , 477 

Fy 

(Lb) 

-3,488 

-3,122 

-3,083 

-3 , 346 

Kxx 


2.352 

2.267 

2.329 

2.344 

E 


-1.461 

-1.378 

-1.280 

-1.397 

iSS 

ESMHH 

-1.998 

-1.874 

-1.871 

-1.961 


(10 6 Lb/in) 

1.573 

1.481 

1.406 

1 .564 

Bxx 

(Lb/ in) 

232.08 

234.79 

269.01 

274.46 

B * y 

(Lb/in) 

-175.53 

-175.87 

-194.38 

-199.65 

Byx 

(Lb/in) 

-174.78 

-174.10 

-192 . 40 

-200 .56 

Byy 

(Lb/ in) 

173.87 

173.79 

187 . 57 

196 .53 
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Figure 38. Comparison to San Andres Five-Pad Bearing at Concentric Position 
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Figure 39. Comparison to San Andres Five-Pad Bearing With 40% Eccentricity) Between Pockets 
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Figure 40. Comparison to San Andres Five-Pad Hearing With 40% Eccentricity Over Pocket 
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At the concentric position, bearing flow requirements calculated by ICYL is 42 versus about 
445/min reported by San Andres. Figure 41 and Figure 42 are plots comparing the direct and 
cross-coupled stiffness coefficients and Figures 43 and 44 compare the direct and cross- 
coupled damping coefficients, respectively, versus eccentricity ratio. In general, ICYL 
predicts about 35% higher direct stiffness, 10% lower cross-coupled stiffness coefficients, and 
15% lower direct damping at the concentric position. With increasing eccentricity ratio, the 
coefficients are observed to behave similarly and some of the discrepancies decrease. The 
cross-coupled damping coefficients with ICYL are equal in magnitude, opposite in sign and 
zero at the concentirc position, as is expected with an incompressible fluid. San Andres non- 
zero concentric value (60 kN-s/m) is due to the fluid compressibility in the pocket. Figure 45 
shows the critical mass versus eccentricity. The concentric value of 119 kg shows better stability 
than predicted by San Andres, which predicts an unstable bearing with a 30 kg mass. 

The analysis of San Andres includes the effect of fluid inertia in the film as well as some 
special effects inside the pocket, such as fluid compressibility and a one-dimensional 
circumferential pressure rise downstream of the orifice. There is also a slight difference in 
friction law used: MTI's analysis follows the formula derived by Nelson [7] for the Moody 
diagram, in which the term containing the Reynolds number is raised to the 1/3 power while 
San Andres uses the same formula with the power changed to 1/2.65 for a more restricted 
range of Reynolds numbers. 

The above comparisons should provide reasonable verification, as the only discrepancies between 
the results can be explained by the different friction and inertia models between the codes. 
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Figure 41. Comparison of Direct Stiffness Coefficient 
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Figure 42. Comparison of Cross Stiffness Coefficients 
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Figure 43. Comparison of Direct Damping Coefficients 
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Figure 44. Comparison of Cross Damping Coefficients 
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Figure 45. Critical Mass versus Eccentricity Ratio 
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